% This function computes the standard errors for the GMM estimates in
% the VAR(1) model when accounting for measurement errors in the factors

function out = getAVarTheta2_threshold_macro_spec(theta2,VarU,...
    Cov10U,Cov01U,x2,macros,econAct,bandwidthT,thresholdLevel,switch_h0,switch_hx)

% Computing the Jacobian matrix
epsValue        = 1D-7;
nx              = size(x2,1)+size(macros,1);
% Remove non-estimated values in theta2
if any(switch_h0 == 0)
    for i=1:nx
        if switch_h0(i,1) == 0
            theta2 = rmfield(theta2,['h0',num2str(i),'_2']);
        end
    end
end
if any(any(switch_hx == 0))
    for i=1:nx
        for j=1:nx
            if switch_hx(i,j) == 0
                theta2 = rmfield(theta2,['hx',num2str(i),num2str(j),'_2']);
            end
        end
    end
end    
namesTheta2 = fieldnames(theta2);
numTheta2   = size(namesTheta2,1);

[~,moments] = objectFuncStep2Threshold_macro_spec(struct2array(theta2)',fieldnames(theta2),VarU,...
    Cov10U,Cov01U,x2,macros,econAct,thresholdLevel,switch_h0,switch_hx);
T = size(moments,2);
%Q = mean(moments,2)'*mean(moments,2)
R = zeros(numTheta2,numTheta2);
for i=1:numTheta2
    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) + epsValue;
    [~,moments_p_eps] = objectFuncStep2Threshold_macro_spec(struct2array(theta2Grad)',fieldnames(theta2Grad),...
        VarU,Cov10U,Cov01U,...
        x2,macros,econAct,thresholdLevel,switch_h0,switch_hx);

    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) - epsValue;
    [~,moments_m_eps] = objectFuncStep2Threshold_macro_spec(struct2array(theta2Grad)',fieldnames(theta2Grad),...
        VarU,Cov10U,Cov01U,...
        x2,macros,econAct,thresholdLevel,switch_h0,switch_hx);
    
    R(i,1:numTheta2) = (mean(moments_p_eps(1:numTheta2,1:T-1),2)...
        -mean(moments_m_eps(1:numTheta2,1:T-1),2))'/(2*epsValue);
end

% The Newey-West estimate of S with lags controlled by bandwidthT (Robust to auto-correlation and
% heteroskedasticity)
GAMMA0 = zeros(numTheta2,numTheta2);
for t=1:T-1
    GAMMA0 = GAMMA0 + moments(1:numTheta2,t)*moments(1:numTheta2,t)'/(T-1);
end
S      = GAMMA0;
GAMMAk = zeros(numTheta2,numTheta2,bandwidthT);
for k=1:bandwidthT
    for t=1+k:T-1
        GAMMAk(:,:,k) = GAMMAk(:,:,k) + moments(1:numTheta2,t)*moments(1:numTheta2,t-k)';
    end    
    GAMMAk(:,:,k) = GAMMAk(:,:,k)/(T-1-k);
    S = S + (1-k/(1+bandwidthT))*(GAMMAk(:,:,k) + GAMMAk(:,:,k)');
end

% The asympotic covariance matrix 
invS             = S\eye(size(S,1)); %inv(S); 
tmp              = R*invS*R';
invTmp           = tmp\eye(size(tmp,1));
%out.VarRobust   = 1/(T-1)*inv(R*invS*R'); 
out.VarRobust    = 1/(T-1)*invTmp; 
out.names        = namesTheta2;
